#include <iostream>
#include <fstream>
#include <math.h>
using namespace std;
#include "network.h"
#include "NODES.h"
#include "SYSTEM.h"
#include "PARA.h"

SYSTEM::SYSTEM() {
    nodes_list=NULL;
    list_len=0;
    }
void SYSTEM::CREATE() {
    if (nodes_list) {
        delete[] nodes_list;
        delete[] nodes_tmp;
        network.~NETWORK();
        para.~PARA();
        SYSTEM();
        }

    para.CREATE();
    if(para.td) { //是否有延迟
        }
    else {
        list_len=2;
        nodes_list=new NODES[list_len];  //存放系统状态,以及计算的中间量2，3，4，5分别为k1,k2,k3,k4
        }
    for (int i=0; i<list_len; i++)  nodes_list[i].CREATE(para.N,para.DOF);
    nodes_list[0].initial(para.tag_initial,para.para_initial);  //初始化系统状态

    nodes_tmp=new NODES[4]; //存放系统临时状态k1,k2,k3,k4
    nodes_tmp[0].CREATE(para.N,para.DOF);
    nodes_tmp[1].CREATE(para.N,para.DOF);
    nodes_tmp[2].CREATE(para.N,para.DOF);
    nodes_tmp[3].CREATE(para.N,para.DOF);

    network.CREATE(para.para_network_tag,para.para_network);
    }



SYSTEM::~SYSTEM() {
    if (nodes_list==NULL) return;
    delete[] nodes_list;
    delete[] nodes_tmp;
    SYSTEM();
    }


void    SYSTEM::function(NODES &y, NODES  &x, double t) { //无延迟的运动方程
    double sumx,sumy,mo;
    for (int i=0; i<para.N; i++) {
        sumx=sumy=0;
        for(int j=1; j<=network.LINK[i][0]; j++) {
            sumx+=(x.STATE[0][network.LINK[i][j]]-x.STATE[0][i])*network.MATRIX[i][network.LINK[i][j]];
            sumy+=(x.STATE[1][network.LINK[i][j]]-x.STATE[1][i])*network.MATRIX[i][network.LINK[i][j]];
            }
        mo=x.STATE[0][i]*x.STATE[0][i]+x.STATE[1][i]*x.STATE[1][i];
        y.STATE[0][i]= para.u*x.STATE[0][i]  -network.FEATURE[i]*x.STATE[1][i] -
                       mo*x.STATE[0][i]  +sumx    +para.f*cos(para.o*t);
        y.STATE[1][i]= para.u*x.STATE[1][i]  +network.FEATURE[i]*x.STATE[0][i] -
                       mo*x.STATE[1][i]  +sumy   +para.f*sin(para.o*t);
//    cout<<y.STATE[0][i]<<"   "<<y.STATE[1][i]<<endl;
//    cout<<"1st: "<<para.u*x.STATE[0][i]<<"  2nd: "<<-network.FEATURE[i]*x.STATE[1][i]<<"  3td: "<<mo*x.STATE[0][i]<<" 4st: "<<sumx<<endl;
//    cout<<"1st: "<<para.u*x.STATE[1][i]<<"  2nd: "<<network.FEATURE[i]*x.STATE[0][i]<<"  3td: "<<mo*x.STATE[1][i]<<" 4st: "<<sumy<<endl;
        }
    }

void   SYSTEM::solve_RungeKutta() { //无延迟龙格库塔
    ofstream file(para.filename().c_str());
    int now=0,future=1,tmp=1;  //指针，代表当前、下一步系统状态以及一个临时交换变量
    double h=para.h;  //时间步长

//  file<<"t  ";
//  for(int i=0; i<para.N; i++)
//    for(int j=0; j<para.DOF; j++)
//      file<<" N"<<i<<"D"<<j<<"  ";  //标题
//  file<<endl;
//
//  file<<para.t0<<"  ";
//  nodes_list[now].print_state(file,'s');

    cout<<para.filename().c_str()<<endl;
    double t=para.t0;
    cout<<"总运行步数："<<para.step<<endl;
    for (int i=1; i<para.step; i++) {
        t=t+h;
//        cout<<i<<endl;
        function(nodes_tmp[0],nodes_list[now],t);

        for(int n=0; n<para.N; n++)
            for(int d=0; d<para.DOF; d++)
                nodes_list[tmp].STATE[d][n]=nodes_list[now].STATE[d][n]+nodes_tmp[0].STATE[d][n]*h/2.0;
        function(nodes_tmp[1],nodes_list[tmp],t+h/2.0);

        for(int n=0; n<para.N; n++)
            for(int d=0; d<para.DOF; d++)
                nodes_list[tmp].STATE[d][n]=nodes_list[now].STATE[d][n]+nodes_tmp[1].STATE[d][n]*h/2.0;
        function(nodes_tmp[2],nodes_list[tmp],t+h/2.0);

        for(int n=0; n<para.N; n++)
            for(int d=0; d<para.DOF; d++)
                nodes_list[tmp].STATE[d][n]=nodes_list[now].STATE[d][n]+nodes_tmp[2].STATE[d][n]*h;
        function(nodes_tmp[3],nodes_list[tmp],t+h);

        for(int n=0; n<para.N; n++)
            for(int d=0; d<para.DOF; d++)
                nodes_list[future].STATE[d][n]=nodes_list[now].STATE[d][n]+(nodes_tmp[0].STATE[d][n]+nodes_tmp[1].STATE[d][n]*2.0+nodes_tmp[2].STATE[d][n]*2.0+nodes_tmp[3].STATE[d][n])*h/6.0;


        if (i>para.step*0&&!(i%5)) { //仅输出一部分数据
            file<<t<<";";
            nodes_list[future].print_state(file,'s');
            }

        tmp=now;
        now=future;
        future=tmp; //总是保证暂时没用的future作为临时变量
        }
    file.close();
    }



//void SYSTEM::run() {
//    network.create(para.para_network_tag,para.para_network);
//string a("net.dat");
//ofstream file(a.c_str());
//network.print_network(file);
//file.close();
//    solve_RungeKutta();
//    }
